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ABSTRACT 


Hydrographic surveys for nautical charting contain many 
discrete data points. Analytical models for ocean bottom 
topography could save computer storage and reduce the com- 
plexity of automating the nautical charting process, but 
they must meet stringent accuracy requirements. Polynomials, 
double Fourier series, finite elements, Duchon's analysis, 
Shepard's formula and Hardy's multiquadric analysis were 
investigated as possible modeling techniques, Multiquadric 
analysis in which the surface is represented by en analyti- 
cal summation of mathematical surfaces such as cones and 
hyperboloids was the only method found to be suitable. ån 
iterative method of model point selection was found to give 
ne best results. Smooth and unambiguous junctions of 
adjacent models were made by using a Hermite polynomial 
weighted sum of overlapping areas, Highly irregular surfaces 
can be represented by about 20% of the original survey data 
points; more regular bottom topography can be represented 


by a smaller percentage. 
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IVI SDG 


A. BACKGROUND 


The ocean bottom is aå continuous but generally irregular 
surface, In the deep oceans there are vast areas of abyssal 
plains interrupted by mid-ocean ridges, sea mounts and con- 
tinents. The continental shelves and coastal areas vary 
from smooth flat bottoms to highly irregular surfaces with 
deeply gouged glacial troughs or coral and rock pinnacles. 
Many geological formations which are found on land such as 
canyons, mountains, domes, faults, etc., are also found on 
the continental shelves. The shape of the ocean bottom is 
difficult to determine since it cannot be seen or photo- 
graphed except in very shallow areas and, direct measurement 
requiring occupation of the ocean bottom is costly and often 
impossible. 

There are many reasons for which the shape of the ocean 
bottom must be known. Historically, safety of navigation 
has been the most urgent reason. Nautical charts are com- 
piled from many sources to aid the navigator. These charts 
depict the coastlines and ocean bottom features using con- 
tour lines and selected depths. 

The primary sources of depth data for nautical charts 
are hydrographic surveys. These surveys represent ocean 


bottom topography by discrete data points which are defined 
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by geographic position and depth below a specified water 
level datum. Until the mid-twentieth century, these depths 
were determined by lowering a weight on a calibrated line 
until it touched bottom, The vessel position was usually 
determined by measurements with sextants. Using these 
manual methods, data acquisition was very slow and only a 
minute percentage of the bottom was sampled. There were 
many sources of error in the observational procedures. A 
typical survey had a few hundred data points from which the 
surface shape between points had to be inferred. Data pro- 
cessing was easily handled by manual methods. More recently, 
electronic positioning equipment and depth sounding instru- 
ments have been used in semi-automated and automated systems. 
These systems allow almost continuous sampling of the ocean 
bottom along the vessel track. They have increased the 
accuracy of the data and the completeness of bottom coverage. 
As a result, depths need to be inferred between vessel tracks 
but not along the tracks. A typical survey of this type 
contains between 2,000 and 20,000 data points. These sys- 
tems increased the data acquisition rate to such an extent 
that manual data processing methods could not xeep up with 
data acquisition. Computer aided systems for processing 
and verifying the data were developed in the 1960's and 
1970's. 

Producing a nautical chart requires compilation of 


many hydrographic surveys, shoreline manuscripts, and other 
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documents. This remained a manual process until the mid 
1970's. At this time, the National Ocean Survey (NOS) of 
the United States National Oceanic and Atmospheric Adminis- 
tration (NOAA) began development of a computer assisted 
chart compilation and production system (Moses and Passauer, 
1979). This system requires on-line storage and manipu- 
lation of large blocks of discrete point data from hydro- 
graphic surveys, The density of these data from modern 
surveys make this a complex and costly process. 

In an effort to produce one hundred percent bottom 
coverage for critical areas, multi-beam sounding systems 
(Hopkins and Mobley, 1978), airborne laser depth measuring 
systems (AVCO Everett Research Laboratory, 1978), and airborne 
water penetrating photography systems (Zeller, 1976) have 
been developed. Some of these systems have proved that one 
hundred percent bottom coverage is feasible. They have 
also created another problem concerning representation of 
the data and its use in the compilation of nautical charts. 
The data from the multi-beam sounding systems for a typical 
survey would be equivalent to several hundred thousand dis- 
crete data points. Data from a laser system would be even 
more dense. The photogrammetric method uses stereographic 
images produced from aerial photographs. This can be con- 
sidered to be truly continuous data, but such data is diffi- 
cult to represent in a digital computer. The usual method 
to represent this data is to select the most representative 


and most critical depths for use as if they were from a 
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conventional survey. For a bottom with little relief, 
this method is satisfactory but as bottom relief increases, 


considerable detail and completeness is lost. 


B. MATHEMATICAL MODELS FOR OCEAN BOTTOM TOPOGRAPHY 


The density of data from modern hydrographic surveys 
has made the automation of chart compilation difficult. A 
possible solution to this problem which is investigated by 
this thesis is the use of a surface defined by an analytical 
expression to approximate the ocean bottom topography. 

Such a mathematical model would be used to compute a depth 
at any geographic position within the bounds of the model. 
In order to be useful, such a model must require consi- 
derably less data storage for the perameters which define 
the model than was required by the original set of discrete 
points. 

The accuracy of the model is of utmost importance, The 
United States government can be held liable for vessel 
groundings or accidents at sea which are due to inaccurate 
charts. Special Publication 44 of the International Hydro- 
graphic Bureau (1968) states the accuracy specifications 
recommended for hydrographic surveys. The depth measure- 


ment specifications are listed in Table I. 
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Table I - Depth Measurement Specifications 
Recommended by the International Hydrographic Bureau 


Depth Allowable error 
0-20 meters (0-11 fathoms) 0.3 meters (1 foot) 
20-100 meters (11-55 fathoms) 1.0 meters (0.5 fathoms) 
Deeper than 100 meters 1% of depth 


The Hydrographic Manual of the National Ocean Survey 
(Umbach, 1976) adds that accuracies attained for all hydro- 
graphic surveys conducted by the National Ocean Survey shall 
equal or exceed the specifications recommended by the Inter- 
national Hydrographic Bureau., These standards do not 
necessarily apply directly to the accuracy requirements 

for a mathematical model of the bottom, but they are good 
reference figures. 

Solution of the dense data problem for nautical charting 
was the primary motivation for the investigation, but there 
are other uses for models appraximating ocean bottom topo- 
graphy. Many coastal processes are closely related to 
bottom topography. These include wave height, wave refrac- 
tion, energy dissipation, wave runup, storm surge and 
beach erosion. Design of offshore structures requires 
input of bottom characteristics. Subsurface, as well as 
surface navigation, could be aided by an ocean bottom model 
stored in an onboard computer. The accuracy requirements 
and model scales for these applications would be different 


but the modeling methods could be the same, 
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C. SCOPE OF WORK 


There are several ways to represent surfaces by mathe- 
matical expressions. Those that seemed most applicable 
to the problem are discussed in Section II. Three of the 
models were chosen for experimental analysis. Portions of 
four hydrographic surveys conducted by the National Ocean 
Survey were used as experimental data sets for this analysis. 
These data sets represent a variation from extreme bottom 
relief to a very flat bottom. The models developed for 
these areas were analyzed quantitatively by comparing 
observed survey depths and computed model depths at the same 
location. Qualitative comparisons of depth contours from 
the two sources were also made. For each type of model, 
the input parameters were varied to investigate minimum 
requirements for a good representation. 

Determining the exact location of the shoreline and 
other boundaries is an important part of any survey, but 
including this in the models is beyond the scope of this 
investigation. All the areas used for experimentation 


were restricted so that they do not include shoreline. 
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II. SURFACE MODELING METHODS 


Analytical expressions have been used previously to 
approximate topographic surfaces, Some techniques used in 
map analysis are also applicable to the problem and there 
are some appealing methods which have been used for other 
surface approximations but not for terrain models. None 
of these methods have been used to represent hydrographic 
surveys. Ocean bottom topography is often similar to land 
topography but the research on terrain models has generally 
been for small scale large area maps. The large scale 
hydrographic surveys which must represent detail on the 
order of tenths of fathoms or feet are quite different than 
those large area maps, so modeling techniques which are good 
for small scale terrain models may not be appropriate for 
hydrographic survey modeling. Some important properties of 
the methods which must be considered aside from accuracy are: 


e ease of computation - Must a large system of equations 


be solved to develop the model? 

e dependence of horizontal scale - Hydrographic surveys 
end marine charts of different scales often overlap or 
are adjacent, For this reason, it is not good if the 
accuracy of a modeling method varies with horizontal 
distance scale, 

e global versus local models - A global model represents 


a large area with a single expression. A local modeling 
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method represents many adjacent small areas with 

many corresponding expressions. Generally, there 

is more computation involved in global methods, whereas, 
local modeling requires more data searching to find 

the appropriate local parameters. Global models 

which attain significant data storage savings are of 
particular interest in this study. 

e interpolation versus approximation - Interpolation 
methods generate a surface which fits some data points 
exactly and is used to interpolate between those points 
for surface values at other positions, Approximation 
methods generate a surface which approximates all the 
data but may not fit any data points exactly. A "best 
fit" by some criteria such as least squares is usually 
used. Approximation methods may not represent the 
least depth in an area accurately or they may move the 
position of peaks and deeps significantly. It is impera- 
tive that the model can be controlled to represent criti- 
cal data points exactly. Interpolation methods are thus 
more appropriate for this application. The data points 
piri aro selected for interpolation will be called 
model points in iNet One center. Quite often they 
are Significant data points such as a least depth or 
an area of slope change. 

The following sections discuss methods and previous 


research which are applicable to the problem. 
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A. POLYNOMIALS 


Czegledy (1977), Hardy (1971), Krumbein (1966), and 
Whitten (1970), discuss the use of polynomials for surface 
representation, A polynomial mapping equation of two inde- 
pendent variables with a specified degree can be produced 
which fits a few data points exactly or approximates all 
the data in a least squares sense. In either case, the sys- 
tem of equations which must be solved becomes ill-conditioned 
as the degree of the polynomial increases. This can De 
alleviated by using orthogonal polynomials. In the method 
of orthogonal polynomials, a collocated series of inde- 
pendent surfaces, linear, quadratic, cubic, etc., is generated. 
The summation of these surfaces is the mapping equation which 
defines the model. Increasing detail is gained by solving 
for and adding the surface of next higher order. This method 
has proven useful for trend analysis of maps. However, it 
has been rejected by some investigators for applications 
requiring more accuracy. The reason as stated by Hardy (1971) 
is that the "ordinary collocated polynomial series is 
unmanageable in representing the sometimes rapid and sharp 
variations of real topographic surfaces." Requiring a high 
degree polynomial to fit closely spaced irregular surface 
paints in one area causes significant invalid variations in 
other areas. To avoid these problems, low degree poly- 
nomials have been used in a local approximation mode with 


success, but this does not produce a global surface model. 
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B. DOUBLE FOURTER SERIES 


The double Fourier series model is discussed by James 
(1966) and Krumbein (1966). It is produced by a series of 
independent harmonic surfaces having wave forms of diminishing 
wave length as the order of the surface increases. This 
technique has proven valuable for trend analysis particularly 
when the surface features show oscillating patterns. Unfor- 
tunately, the models require high order surfaces to repre- 
sent sharp terrain features. Such surfaces produce oscillations 
with large variations between data points and have many of 


the same drawbacks as the collocated polynomial series, 


C. FINITE ELEMENTS 


Gold, Charters and Ramsden (1976) discuss a method of 
surface representation in which a system of triangles with 
data points at the vertices is imposed on a surface. An 
interpolating function is used to estimate the surface in 
each triangular element. The interpolant is developed so 
that the surface passes through the vertices and makes a 
smooth transition from one triangle to the next. 

Peucker, Fowler, Little and Mark (1977) have developed 
a similar system of surface representation by Triangulated 
Irregular Networks (TIN). Rather than a smooth interpolant, 
the TIN system uses the planes defined by the three function 
values at the vertices of each triangle to represent the 


surface. Considerable work has been done on automated 
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techniques for selecting appropriate points to be used 

for vertices and on development of data structures for 
storage of the vertices, neighboring points, and neigh- 
boring triangles. The TIN system was developed specifically 
for digital representation of topographic surfaces. 

Finite element systems such as these are local methods. 
Detail can be easily incorporated into the model by adding 
points where required without affecting the model elsewhere, 
Very little computation is required but searching the data 
structure to find the appropriate element is necessary. 

Such systems are generally independent of scale unless a 
scale dependent interpolant is used. A single expression 
which represents the surface is not generated by these 


methods. 


D. SHEPARD'S FORMULA 


Shepard's method as described by Poeppelmeier (1975), 
Barnhill (1977) and Franke (1979), nas been widely used to 
interpolate random data but has never been used for topo- 
graphic surface representation. The model is produced 
by taking a weighted average of the model points to inter- 
polate the surface value at other points, 


Shepard's formula is expressed by 


2. e e 
i wifi if d; + O for all i 
2. 
a i “i 
Í; LE d; = O for any i 
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where the Í; are the depths at the model points; d; is the 
distance from the ith model point to the point of computa- 
tion; and the weight assigned to each model point, Wis is 
a unction of =. Two such weighting functions used in 
this project R simply the inverse distance (1/d, ) and 
the inverse distance squared AO. 

In this method, all model points contribute to the 
value of f, but the effect of any model point on the inter- 
polant decreases as the distance from that point increases. 
Another appealing feature of this method is that the value 
of f will always be between the minimum and maximum values 
of the model points. 

Franke and Little's modification to Shepard's method 
restricts the weighted summation to only those model points 
within a radius R of the computation point. With this modifi- 
cation, the weighting function approaches zero as the dis- 
tance approaches R and remains zero at distances greater than 


R. The modified Shepard's formula is expressed by 


2 (R-d.), 
en 


Rd; fi 
Z i/+ 
f = oS Ti (2) 
Í; es di = O Lor any i 
or 2 
2 (Reds) £ 
i 2%4 2 i 
i if d; 4 O for all i 
— 
>: (Rudy), 
i 
f = R d; (3) 
fi If d; = O for any 1 
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where 
R-d Erla 


(R-d,), = (4) 
pr O ford, >R 


The weighting functions 1/4, and (R-d;), produce surfaces 
š i 
with cusps at the model points. The weighting functions 


2 ' 
1/47 and (R-d;), produce surfaces with flat spots at those 


R à 


points. For higher order functions of l/d; these flat spots 


increase in size and the slopes between them become steeper. 


These properties are shown in two dimensions in Figure l. 
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Figure 1 - Shepard's Formula with Various 
Weighting Functions 
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These formulas do not require solution of systems of equa- 
tions and are easily modified by simply adding significant 
data points without recomputing any coefficients. They are 
independent of scaling, global in nature, and the computa- 


tion is very simple. 


E, DUCHON'S METHOD 


The method of Duchon (1976) which was developed as thin 
plate surface theory is described by Meinguet (1979) and 
Harder and Desmarais (1972). It has never been used for 
topographic surfaces but has been used for other surface 
analyses. To develop this model, individual surfaces called 
basis or kernel functions, which are centered at the model 
points, are summed to yield a global surface. There is a 
coefficient associated with each kernel function which 
determines the magnitude of the effect of that kernel func- 
tion on the total surface. 

The expression for the model is 

n 
f = Z Š: Fa, Y, Xi J + Ay + boxe eds Y (5) 
where n is the number of basis functions and model points 
used, The last three terms represent a plane which is also 
added into the model. The n+3 coefficients Cas dll; e wis 
Az y A, and Az are determined by solving the following system 


of n+3 equations. 
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where Íj Los —--, f, are the surface values at the model 
points (XY) (X, Yo), Zu, (X vor 


Duchon used two basis functions 


Hi Xis Ti) si sé 
and 2 Ci) 
Peet. Xis Yi) — Ma log dy 
where 2 2, 4 (8) 
d; =((X-X;)* + (Y = ¥,)°) $, 


d; is the horizontal distance from each model point to the 
point of computation, 

Duchon's method using the above basis functions is 
independent of scale. During experimentation, a third basis 
function 


uy s Opa 


) = d; 108 d; (9) 
was also used. The models using this latter basis function 


are dependent on scale. 
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F, HARDY'S MULTIQUADRIC ANALYSIS 


Multiquadric analysis, as discussed by Hardy (1971, 
1972a, 1972b, 1975 and 1977) resulted from a search for a 
satisfactory and efficient method to represent topography 
by an analytical model. As suggested by its name, the 
method consists of summing many quadric surfaces (cones, 
hyperboloids, paraboloids, etc.), each associated with a 
model point, to obtain a global surface. Superficially, 
this method is similar to Duchon's method except that the 
kernel functions are quadric surfaces and the additioanl 
three terms are not used. The expression for this model 
LS 


f = £ > en E eo Y,)] (10) 


5 = 


where f is the surface value at the point (X, Y); Q is the 
quadric surface or kernel function; (X; 9 Y), i=l, ---, n 
are the model points at which the kernels are centered; 
and Cis i=l, ---, n are coefficients assigned to each sur- 
face. 

The following system of equations is used to solve for 


the unknown coefficients. 


n 
i E (Xp, Ty Xi; J 


pa 
l 


(135) 
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Theresulting surface will fit the data exactly at the model 


points (X; 9 Y f.), i=l, ---, n. 


317 
Hardy (1971) found that the quadrics which yield the 
best reults are hyperboloids, cones and inverse hyperboloids. 

A hyperboloid is represented by 
Sr r +) 2) 
GS i i i 
Cones are special cases of hyperboloids where Ó 18 zero: 


OX, Yr Xp i 


1 
d, is the distance from the point of computation (X OSEO 


the center point of the each quadric surface (X, Yi). 


Inverse hyperboloid kernels are expressed by: 


(X, Y, X, Y) = ((X-2,)? + (7-1,)? + 82%, (14) 


The magnitude of the coefficient Ci determines the steepness 
of the cone or hyperboloid. The sign of C; determines 
whether the surface is oriented upward or downwerd. The 
magnitude of Ô determines how flat the hyperboloid is at 


its center. These properties are shown in Figure 2, 
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Figure 2 = Hyperboloid Kernels 
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For an inverse hyperboloid ð determines the peakedness of 
the surface at its center point (Figure 3). C; simply 
represents a multiplicative constant which scales the size 
of the surface and specifies its orientation upward or 


downward. 





Figure 3 - Inverse Hyperboloid Kernels 


For all Ci and Ò the inverse hyperboloid approaches zero 
as the distance increases, There is an inflection point 
at d = N=. 
The way several quadric surfaces sum to form the global 


surface can be seen in two dimensions in Figure 4. 
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Figure 4 - Quadric Summations (hyperboloid kernels) 


Multiquadric modeling is independent of horizontal 
scaling so long as Dis linearly related to the horizontal 


distance scales. 
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III. HYDROGRAPHIC SURVEYS 


The data used in this project were from hydrographic 
surveys conducted in the late 1970's by the National Ocean 
Survey. The methods and procedures used were typical auto- 
mated survey procedures as documented by Umbach (1976) and 
Wallace (1971). A brief description of these procedures 
followed by more specific information on each data set 


follows, 


A. GENERAL DESCRIPTION 


Safety of navigation is the primary purpose for which 
hydrographic surveys are accomplished, The data is acquired 
by running sounding vessels in parallel or nearly parallel 
tracklines on the ocean surface and taking depth measure- 
ments along these lines at evely spaced intervals. Cross- 
lines are run at large angles to the main system of lines 
as a gross check on the validity of the data. When indica- 
tions of critical bottom detail are found, development lines 
are run at closer spacing to determine the least depth and 
verify the nature of the feature. The depths are plotted 
on a survey sheet, a sample of which is shown in Figure 5, 

There are many properties of a survey which affect its 
usefulness. Those most important to this project are survey 


scale, horizontal positioning accuracy and depth accuracy. 
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Figure 5 - Portion of a Hydrographic Survey Sneet 
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in Survey Scale 


The survey scale is the ratio of distance on the 
survey sheet to the corresponding distance on the earth. 
The scale chosen for a survey depends on "the area to be 
covered and the amount of detail necessary to depict ade- 
quately the bottom topography and portray the least depths 
over critical features." (Umbach, 1976) The survey scale 
is usually at least twice as large as the scale of any 
chart published for the area. Large scale surveys cover 
less area than small scale surveys but greater detail can 
be represented. For this reason, large scale surveys are 
conducted in harbors, anchorages, restricted navigable 
waterways, and areas where dangers to navigation are 
numerous. Areas with considerable detail are the most 
difficult to adequately represent by a mathematical expres- 
Sion. Three of the four data sets used in this project 


were from large scale surveys. 


2. Horizontal Position Accuracy 
Umbach (1976) specifies that plotted positions, 


"whether observed by visual or electronic methods, combined 
with plotting error shall seldom exceed 1.5 mm (0.05 in.) 
at the scale of the survey," On a 1:5000 scale survey, 

the position of each sounding should thus be represented 

to within 7.5 meters of its actual position on the earth, 
This is important in evaluating a mathematical model. One 


of the data sets had some very steep slopes, where an error 
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of a few meters in positioning would produce a depth 
variation of several fathoms. For areas such as these, a 
much greater depth discrepancy between the model and the 


survey data should be tolerable. 


3. Depth Accuracy 


As seen in Figure 6, there are many components 
that make up the depths represented in a hydrographic 
survey. In addition to the depth recorded by the sounding 
instrument, there are corrections for velocity of sound in 
the water column, the stage of the tide, and the dynamic 
vessel draft. Sometimes surveys have slight inconsistencies 
where data from two different vessels or two different days 
are adjacent or intermixed. These might be due to changes 
in the water column structure that affect the velocity of 
sound, an error in determining offshore tide corrections 
from tide gages near the shore, unrecorded changes in vessel 
speed affecting the dynamic vessel draft, a slight systematic 
error in vessel positioning, etc. Even more critical is 
the effect of waves on the sounding vessel. Small vessels 
tence vertical position rapidly as waves pass while the 
instruments record the depth of the water column below the 
vessel. This depth is too great if the vessel is on a 
wave crest and too small if the vessel is in a trough. The 
angular orientation of the vessel is also affected by waves. 
If the vessel rolls to an angle greater than the sounding 


beam width, the depth recorded may not be under the vessel 
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Figure 6 - Components of a Depth Measurement 
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bu off to the side. In order to correct for this, the 
echograms were manually scanned, wave action was visually 
meaned out of the record, and depths that were automatically 
acquired with an error greater than the recording interval 
were rerecorded. On days of moderate to heavy wave action, 
this procedure leads to an inordinate amount of manual 

work introduced into an otherwise automated system. Table II 
indicates the depth recording and correction intervals used 
by NOS, Note that in many cases soundings from 0-20 fathoms 
need only be recorded to the nearest whole foot or nearest 
half of a fathon. 

Althoush depth measurement errors can exist on all 
surveys, they are more apparent in areas of flat regular 
bottom. If two adjacent soundings each have nearly a foot 
of error of opposite sign, this will appear as a sharp dis- 
continuity on a flat bottom whereas it will hardly be 
noticed on a steep slope. For this reason, models for 
areas of flat regular terrain when compared with the survey 
data may show some relatively large differences due to the 


data acquisition procedures. 


B, DATA SETS 


Four deta sets were used for analysis in this project. 
They were specifically selected for the variety of bottom 
topography which they represent, 
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i Monterey Bay, California 


This data set was taken from survey registry number 
H-9808. It was conducted in 1979 by NOAA Ship DAVIDSON and 
Naval Postgraduate School personnel and equipment. It 
covers the southernmost part of Monterey Bay including 
Monterey Harbor. The survey was conducted at a scale of 
1:5000. Only one vessel was used on the portion of the 
survey chosen for analysis. The sounding units are fathoms 
and depths range from 0 to 16 fathoms. The bottom has a 
large amount of detail. It slopes moderately downward from 
the shore and consists generally of mud and sand. In the 
middle there is an area thick with kelp which is attached 
to a rocky irregular bottom. There are a few rocky areas 
in the deeper part as well. Figure 7 shows the bottom con- 
tours in one fathom increments. The scale of the plot has 


been reduced for presentation nerein. 


2. Morro Bay, California 


This data set was taken from survey registry number 
859757. It was conducted in 1978 by the NOAA Ship FAIRWEATHER, 
It covers a small part of Morro Bay and some navigable water- 
ways open to the bay. The survey was conducted at a scale 
of 1:5000. The sounding units are feet and depths range from 
16 to 82 feet in the portion used for analysis. Figure 8 
(reduced scale) anes the bottom contours in three foot incre- 
ments. There is one major feature near the center and con- 
Siderable irregularity in the northeast corner of the area. 


Otherwise, the bottom slopes gently offshore. 
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3. Auke Bay, Alaska 


This data set came from a thesis project by Seidel 
(1979), a student at the Naval Postgraduate School, which 
investigated the affects of using multiple sounding beam 
widths for hydrographic surveys. The procedures were some- 
what non-standard since sounding lines were run much closer 
than normal in an attempt to gain 100% bottom coverage. 
Specifications for 1:5000 scale surveys were used but due 
A the dense sounding spacing, it was plotted at a scale 
of 1:2500. The data was incorporated into survey registry 
number H-9818. It was conducted in 1979 by Seidel and the 
NOAA Ship RANIER. It covers a small portion of Auke Bay in 
southeast Alaska. The sounding units are fathoms and depths 
range from O to 24 fathoms. The bottom is mostly mud and 
rock and shows a tremendous amount of variation due to 
glacial action. Very steep slopes are encountered in the 
area. At one point, the depth changes from 7 to 22 fathoms 
in a horizontal position change of only 30 meters. Figure 9 
(reduced scale) shows the bottom contours of the central 


part of the data set in one fathom increments. 


4. Gulf Coast 


The fourth data set was taken from survey registry 
number H-9785. It was conducted in 1978 by the NOAA Ship 
MT MITCHELL at a scale of 1:20000 and covers an area in the 
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Gulf of Mexico off the coast of Louisiana. The sounding 
units are feet and depths range from 29 to 37 feet in the 
portion used for analysis. Figure 10 (reduced scale) shows 
the bottom contours in one foot increments, The bottom is 
generally flat with a very gentle slope, It consists 

mostly of mud and shell fragments. Some of the irregulari- 
ties seen in the bottom contours are in areaswhere the work 
of two vessels overlapped because of crosslines or junctions. 
The flat bottom and small contour increment make these 
irregularities stand out. The survey party reported that 


wave action was also a considerable problem during the con- 


duct of this survey. 
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Figure 10 - Gulf Coast Data Set Contours (feet) 


46 





IV, RESEARCH PROCEDURES 


A, COMPUTER SYSTEM 


All the computer work for this project was done on 
the IBM 360/67 system at the Naval Postgraduate School's 
W. R. Church Computer Center. All programs were written 
in FORTRAN IV. The VERSATEC-07 electrostatic plotter was 
used for all the data and contour plotting, Both IMSL 
(International Mathematical and Statistical Library) 
routines and other library routines were used in the pro- 


grams, 


B. DATA SET PREPARATION 


l. Original Data Condition 
A11 four data sets were supplied by the National 


Ocean Survey on non-labeled unblocked magnetic tapes in the 
NOS standard record format. Positions of all soundings were 
given in terms of latitude and longitude. Corrected soundings 
were supplied to the nearest tenth of feet or fathoms. Each 
data point had a record sequence number assigned, The NOS 
format also included original observed data and all correc- 
tions to it as well as descriptive cartographic codes and 


other information. 
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2. Program TAPCNV - Tape Conversion 


Only corrected position, corrected depth and 
record number identification were required for this pro- 
ject. The program TAPCNV was written to read this data 
from the non-labeled NOS tapes. The geodetic positions 
- were converted to an X-Y plane coordinate system based on 
the Modified Transverse Mercator (MTM) projection (Wallace, 
1971). Double precision computations were used for this 
conversion. The MTM projection gives the positions in 
terms of meters of northing and meters of easting from 
a local origin. This X, Y position was then converted to 
plotter coordinates in terms of inches from the plotter 
‘origin. The record number, depth, geodetic position, 
MIM coordinates and plotter coordinates were blocked and 
recorded on disk and on an NPS tape with standard system 


labels. 


3. Program DATPLT - Data Plotting 


This program was written to display the discrete 
point data on a plotted sheet, Latitude and longitude 
grid intersections at specified intervals were converted 
to plotter coordinates and straight lines were drawn con- 
necting these points to provide the geodetic position 
reference system. Two sheets were plotted with this 
reference grid. On one sheet depths were plotted to the 
nearest tenth. (NOS plots tenths only in shallow water 


when the depth units are fathoms.) Record numbers were 
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plotted on the second sheet. Overlaying the two sheets 
facilitated reference to any particular data RS DAT 
was used to plot the entire surveys as a first step. When 
portions of the surveys were selected for analysis DATPLT 
was used again to select and plot only the data within 


the specified area, 


4. Program CONDAT - Data Contouring 


Part of the data analysis consisted of comparing 
contours of the original data with contours from the model. 
Initially, contours of the survey data were hand drawn - a 
procedure that is somewhat subjective. In order to remove 
as much subjectivity as possible from the analysis hand 
contouring was replaced by machine contouring. The library 
routine CONISD, for contouring irregularly spaced data, was 
used in the program CONDAT, This routine first generates 
triangles with data points at the vertices, By linearly 
interpolating along the triangle sides for the contour 
values, points on each contour are found and connected to 
generate the contour lines, The contours generated in this 
manner were not smooth as would be desirable, but the data 


points were dense enough so that this was not a problem. 


C. MODEL DEVELOPMENT AND ANALYSIS 


Program MODEL was written to do the model development, 


the quantitative analysis, and to aid in the qualitative 
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analysis. To change from one modeling technique to another 
the only modification necessary was the replacement of one 
module, That module contained the routine to develop the 
model by a given method and to compute the depth at any 
point. All the modeling techniques required the selection 
and use of model points from the survey data, The model 
points were specified by record numbers on punched card 
input and the survey data was read from disk and stored in 
memory. The model points were stored in arrays for model 


development. 


J, Coefficient Computation - Subroutine LEQ2S 


The methods of Hardy and Duchon require solution 
of symmetric systems of linear equations to determine the 
model coefficients. The double precision version of the 
IMSL routine LEQ2S was used for this purpose. This routine 
uses symmetric decomposition with iterative improvement 
to solve the systems. Systems of up to 226 equations in 
226 unknowns were solved during the course of this project. 
The model coefficients and respective model points were 


output for analysis. 


2. Quantitative Analysis = Subroutine STAT 


Subroutine STAT was developed to provide a quanti- 
tative analysis of each model, Bach survey depth was com- 
pared with the depth computed from the model at the same 
X, Y position. The root mean square difference, maximum 


positive difference and maximum negative difference were 
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tabulated for each run. A positive difference signified 
that the model depth was deeper than the survey depth; a 
negative difference signified that the survey depth was 

deeper. The root mean square difference is given by the 


expression 
N 2 
MD-SD) È 
à 


RMS difference = A eee (15) 
N 


where MD is the model depth, SD is the survey depth, and 
N is the number of points used for the comparison. Those 
points with a difference greater than 1.5 times the RMS 


difference were listed for menual inspection end anelysis, 


3. Qualitative Analysis - Subroutines SETCON and CONTUR 





The qualitative analysis was accomplished by compari- 
son of model contours with those from the original data. 
In order to produce the model contours subroutine SETCON 
developed e quarter inch grid over the modeled area at the 
scale of the survey. The model depth was computed at each 
grid intersection. These depths and positions were passed 
to the library routine CONTUR for contouring. CONTUR is 
similar to CONISD except that it was written specifically 


for gridded data and runs considerably faster than CONISD. 
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V. RESEARCH RESULTS 


A. SELECTION OF METHODS FOR EXPERIMENTATION 


Of the modeling techniques discussed in Section II, 
Duchon's method, Shepard's formula and Hardy's multiquadric 
analysis were selected for testing. Neither Duchon's method 
nor Shepard's formula had previously been used for terrain 
surfaces, Multiquadric analysis had been used with good 
results for topographic data but not at the scales and 
accuracy requirements necessary for hydrographic survey 
representation. 

The methods of polynomials and double Fourier series 
were not tested. They had proved to be useful for some 
applications such as trend analysis and representation of 
repeated features. The fact that forcing polynomials or 
double Fourier series to fit irregular data in small areas 
produces unwanted irregularity in other areas would seem 
to preclude these methods from producing good results in 
this application. Some methods reduce this effect by using 
local expressions which fit only small areas at a time, but 
this defeats the purpose of generating a global model to 
represent large areas of the data. 

Finite element methods were not tested. They are strictly 
local methods which, in addition to storage of model points, 
require storage of pointers to the neighboring points and 


neighboring triangles of each model point. 
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B. METHOD COMPARISON PROCEDURES 


The Monterey Bay and Morro Bay data sets were used for 
comparison of the methods. The Monterey Bay data set was 
used in the qualitative manner only. The Morro Bay data 
set was used for both qualitative and quantitative compari- 
son. The procedures described below were used for all methods 
in order to make controlled comparisons. 

As a first step, 42 data points on a 6x7 grid were chosen 
from the Monterey Bay data at regular spacing without regard 
to bottom detail. After the models were generated with these 
points, an additional 18 model points were selected in areas 
where more detail was required. Models were then generated 
using the 60 points. The third step was to choose 30 more 
points around the outside of the original area at the same 
Spacing as the original 42 points extending the grid to 8x9, 
Models were generated with the 90 points to determine the 
affect of extending the model area. £ 

Thirty-seven model points were chosen at regular spacing 
from the Morro Bay data set in the first step with that 
data. Thirty and thirty-one additional points were selected 
in the second and third steps for totals of 67 and 98 model 
points. The additional points were all within the original 
area in places where additional detail or accuracy was needed, 
The effect of increasing the model point density was examined 
in this way. The statistical results, as well as contours 


of the Morro Bay tests, were compared. 
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C. RESULTS OF DUCHON'S METHOD 


Equation 5 shows that Duchon's model is produced by the 
summation of a plane and a series of basis functions centered 
at the model points. Duchon's method with the basis functions 


2 


a? and d” logd and a similar method with basis function d logd 


were tested. 


1. General Findings 
For all three basis functions the 42 Monterey Bay 
points produced similar models which showed the general 
trend of the bottom but very little detail. Using the 60 


model points, the basis functions a? and a“ 


logd gave more 
detail and a fair representation of the bottom trends. 
See Figure 11. 

The basis function d logd gave a very poor repre- 
sentation of the bottom when the 60 points were used. In 
an effort to resolve detail in some areas, several points 
were chosen very near each other. This caused steep slopes 
to be generated which extended into areas where there were 
no model points and created invalid peaks and deeps. 

The quantitative results of Duchon's method using 
the Morro Bay data set are given in Table III. Note that 
the model using the basis function d logd became better in 
both RMS difference and maximum difference as the number 
of points was increased. The results didn't change much 
using the basis function a“ logd and they became worse for 


the basis function a2, The contours reflected these 
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Figure 11 - 60 Point Duchon Model of Monterey Bay 
(basis function d2 log d) 
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MA tn the cases of the ad’ and af 


logd basis func- 
tions, the detail was increased in the areas where points 
were added but the statistical results did not improve due 
to greater error in other areas, As shown in Figure 12, 
the model using 98 points with the basis function d logd 
compares very favorably with the original data contours in 
Figure 8. 

TABLE III - Results of Duchon's Method - Morro Bay 


WE Of Maximum Maximum IT of 
Basis model positive negative data 

Function points difference difference difference points 
a? 37 1.20 3.89 Zoe 936 
d“logd 37 da 3,78 6,29 936 
d logd ER 1715 4.32 -6,51 936 
a? 67 I; 6.74 -7.99 936 
d“logd 67 ala 4.09 -4,42 936 
d logd 67 CENT Set -3,15 936 
að 98 2.13 7.75 -17.83 936 
d“logd 98 a 4.19 -6,15 936 
d logd 98 9767 23825 -2.48 9356 

2. Devendence on Scale 


The results of the previous section lead to a ques- 
tion concerning the ability of the basis function d logd to 
produce a much better model than other basis functions in 
one case but not in the other. The reason for this turned 


out to be the scale of the data. The first two basis 
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functions produce models which are independent of scale 
whereas, the third produces models which are not, Initially, 
the Morro Bay data was used at the earth's natural scale in 
meters. At this scale, the reale equations was ill- 
conditioned to such a degree that it couldn't be solved. 
The data presented in the previous section was acquired 
with the horizontal position data scaled to a distance of 
unity on a diagonal from one corner of the area to the 
opposite corner. In this case, the results were good. 
The Monterey Bay data was scaled for a diagonel distance 
of 50. The results in this case were poor. Table IV gives 
the results of some tests run at various scales using the 


Morro Bay data and the basis functions d logd and a“ 


logd. 
The set of 98 model points and 956 data comparison points 


were used for these tests. 
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TABLE IV - Effects of Scale on Duchon's Method - 
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Morro Bay 
Basis function: d logd 
Maximum Maximum 
Diagonal RMS positive negative 
distance difference difference difference 
9. 001 036716 Zelo -2,365 
5,01 0:672 25121 nd 
0,1 0.668 24151 -2,406 
1,0 02071 2,216 =2, 472 
10.0 555827 155966 -20,285 
100,0 17.294 6 208 -198.561 
1000.0 O 7151 zo -2,461 
Basis function: a“loga 
0,001 l 45. 47156 -6.149 
a. Ol 1.157 des -6.143 
0.1 15159 4.193 -6.155 
1.0 1,153 Aue) -6.151 
10,0 1,198 4.189 -6,151 
100,0 15153 4.189 -6.151 
1000.0 1.152 45190 =6,155 





De RESULTS OF SHEPARD'S METHOD 


As defined in equation 1, Shepard's formula is a 
weighted summation of the model point depths. The weight 
assigned to each model point is a function of the inverse 

distance from the point of computation to the model point. 
The weighting functions used in this analysis were: 


a 
d. 


1 


o —+ 
a. È 


i 


2 
E F 
R 4.2 
i 
where (R-d, ), is defined in equation 4. 
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In the modified Shevard's method, a radius of in- 
fluence, R, is used. Rather than choosing this radius 
arbitrarily, a method was used which related R to the den- 
sity of the model points independent of the scale of the 
data. This also allowed variation of R according to the 
average number of model points which would fall within the 
radius of influence. 


The following expression was used for this purpose: 


H CAL" NPPR (16) 
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where D is the maximum distance between any two model 
points, N is the total number of model points, and NPPR 
is the average number of points which should fall within 
the radius of influence. FY is related to the size 

of the area which contains the data points. Dividing 

this by N gives a measure of the average area which could 
be assigned to each point. Multiplying by NPPR gives the 
area which could be associated with that many data points. 
Taking the square root of this gives a radius which would 
define that amount of area centered at the point of compu- 
tation. On the average there should be NPPR model points 
within a distance R from any point of evaluation. The 
tabulated statistical results express the radius in terms 


of NPPR instead of a. 


Ze Inverse Distance Weighting Function 
Table V gives the statistical results of the tests 


using the inverse distance weighting functions. The table 
shows that use of the modified Shepard method improved the 
results considerably. In all cases, the best results were 
obtained by including an average of six model points in 
the radius of influence. The table also indicates that no 
statistical improvement was made by increasing from 37 to 
98 model points. 

The contours produced by this method (Figure 13) 
for both data sets were poor. The basic trend of the bottom 


can hardly be seen. The contours are quite wavy where they 
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TABLE V - Shepard's Formula with Inverse 
Distance Weighting Function - Morro Bay 


Number Maximum Maximum Number 
of model RMS positive negative of data 
points NERI di-feerenceemMid-ierence difference points 
D :7 6 2.06 8.17 Her 936 
Di J 2429 Ð -9,68 936 
Sal All* 10,88 24,43 -25.96 9356 
67 4 2.45 8.50 -7.97 936 
67 2,17 30 -6.94 956 
67 3 20 6,86 -7.61 956 
67 25 3,53 9,61 -12.68 936 
67 A11* 11541 Pe -238.58 936 
98 4 za Salt -3,95 958 
98 6 29.17 Tei -6.73 936 
98 2 222 3.64 7,51 996 
98 25 3.46 126 -12.51 936 
98 56 2.05 15359 -16,50 956 
98 All* 12/20 21399 -32.22 936 


* Shepard's formula =- All model points contributed to 
the weighted average. 
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should be straight. In some cases, peaks or deeps are pro- 
duced at the positions of model points which aren't found 


in the original data. 


3. Inverse Distance Squared Weighting Function 


Table VI gives the statistical results of similar 
tests using the inverse distance squared weighting function. 
There is considerably less variability as NPPR is changed 
using this weighting function. The results are better for 
large NPPR and for the unmodified version, but the best 


results at smaller NPPR did not improve. 


E. RESULTS OF HARDY'S MULTIQUADRIC ANALYSIS 


As indicated in equation 10, Hardy's multiquadric model 
is generated by summing quadric kernel surfaces, each of 
which are centered at model points. Hyperboloids, cones 


and inverse hyperboloids were the kernel surfaces tested. 


l. Determination of Ó 
Both hyperboloids and inverse hyperboloids require 
the parameter Ô (Section II.F). Variation of Ó makes 
considerable difference in the results. The effect of 
any value of Ò on the shape of the quadric surfaces with 
respect to the entire model is related to the scale of 
the model. Hardy (1977) has indicated that the optimum 


value of Ô in his investigations was also related to the 


distance between model points. The following expression 
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TABLE VI - Shepard's Formula with Inverse 
Distance Squared Weighting Function - 
Morro Bay 


e 


Number Maximum Maximum Number 

of model RMS positive negative of data 

‘points NPPR difference difference difference points 
SÙ 4 2805 9.51 -3,90 936 
DI 2.37 8.75 -9.09 996 
on 25 2.36 8.38 -9.57 936 
57 All* 5500 14562 -17.06 936 
67 4 Pa ol 9.14 -8.87 936 
67 J 2:92 6.01 -6,88 956 
67 ca 2.54 c 71 -9,00 956 
67 ALIS 2.24 15.42 -19.90 936 
98 4 2395 350) -9.46 936 
98 6 2.95 023 -7.84 936 
e J 2515 (eus -7.05 956 
98 25 2,206 9596 -3.40 956 
98 56 re 10:97 -11.66 936 
98 All* 6.58 20690 — 4. 956 


* Shepard's formula - All model points contributed to 
the weighted average. 
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was used to relate ð to the average density of the model 


Points: 


6 = è 2- x 0.1 x NPPR (17) 
With this expression, the effect of Ò on models of 
different scales will be similar as long as NPPR de 
same. The tables in the following sections are expressed 
in terms of NPPR instead of the absolute value of 6. 

A cone is a special case of hyperboloid where Ò is 

Zero. The tables for hyperboloid kernels include cones by 
listing NPPR as zero. A zero value of Ó is not valid 


for inverse hyperboloids since the peak of an inverse 


hyperboloid increases to infinity as ð approaches zero. 


2. Inverse Hyperboloid Kernels 
For both data sets when ð was small (NPPR=5), the 


contours showed holes at the model points which were not 
indicated in the original data. The representation of the 
actual surface was very poor. Increasing NPPR to 10, 15 
and 20 gave somewhat better results and the bottom trends 
were evident but the representation was still not good. 
With NPPR greater than 20, very steep slopes were created 
in large areas where no model points were chosen, 

The statistical results using the Morro Bay data 
set are given in Table VII. The results became worse as 
more model points were added, The best results were con- 
siderably poorer than the best results from other methods, 


particularly the maximum differences, 
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TABLE VII - Multiquadric Analysis with 
Inverse Hyperboloid Kernels - Morro Bay 


Number Maximum Maximum Number 

of model RMS positive negative of data 

points rr dimrerence y daiigerenee difference points 
37 5 4.99 4.71 35.48 936 
Dali 10 2,05 eee -21.09 956 
Du 15 1259 4.56 -14.73 956 
57 20 1.43 4.24 -10.90 356 
67 Ð oli 90 -32.06 956 
67 10 2.54 T62 lei 936 
67 15 2.69 17.74 =15., 1 930 
67 20 4.09 20.16 -15.60 936 
98 5 251099 CAD -60.04 956 
98 10 Oto F. Ja -23.96 356 
98 15 2:75 21397 -20.13 936 
Je 20 Saa SO =) 22 256 


ST 





3. Hyperboloid and Conic Kernels 


Hyperbolic and conic kernels were evaluated on 
both data sets and NPPR was varied from zero to 25 for 
each set of model points. With 42 regularly spaced model 
points on’ the Monterey Bay data set, not much detail was 
evident but the general bottom trends were well represented 
for all values of NPPR. Increasing to 60 model points 
gave more variation with NPPR. For NPPx set to zero and 
one the results were very good. See Figure 14. The 
detail was improved and the bottom trends were still accu- 
rate. For NPPR set to 10 and 20, the results became pro- 
gressively worse. Very steep slopes were generated which 
created invalid peaks and deeps in areas where no model 
points were chosen. The reason for these slopes is apparent 
when examining the magnitude of the coefficients, For 
NPPR=1, the mean coefficient magnitude was 0.33; for NPPR=20, 
the mean coefficient magnitude was 151.63. hen NPPR was 
increased to 25, the system of equations became so ill- 
conditioned that it could not be solved. This is due to 
the increased flatness of the hyperboloids when NPPR becomes 
large. In areas where several model points are very close 
in order to represent sharp irregularities in the bottom, 
the flat hyperboloids centered at those points can't produce 
the detail required. Increasing to 90 model points produced 
Similar results. 

The statistical results from the Morro Bay data set 


are given in Table VIII. For small NPPR, the results became 
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Figure 14 - 60 Point Multiquadric Model of Monterey Bay 
(hyperboloid kernels) 
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TABLE VIII - Multiquadric Analysis with 
Hyperboloid and Conic Kernels - Morro Bay 


Number Maximum Maximum Number 

of model RMS positive negative of data 

points NPPR difference difference difference points 
317 O 1915 SIA -6.59 256 
37 1 LA 6.00 -6.76 936 
DI 10 j Se -6.92 936 
DI 25 1.24 3.90 -6.17 936 
67 O 0.75 Sol -3.09 950 
67 1 0278 3.99 -3.01 936 
67 10 2.10 15.41 -6.79 936 
67 25 12504 62590 -44.92 936 
98 O CCS 1,92 -2,14 936 
98 I 0.68 2.06 -2,12 936 
98 10 2,54 12512 -14.11 956 
98 20% 14,44 115796 -43,61 9250 


* system couldn't be solved for NPPR=25 
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continually better as more points were added for greater 
detail. As the model points became more dense, the best 
results were acquired by using conic kernels (NPPR=0). 
For the original 37 regularly spaced model points, the 
best results were for small but non-zero NPPR. The con- 
tour comparisons reflected the model quality demonstrated 


by the statistics. 


F, SUMMARY 


A graphical comparison of the statistical results of 
the methods is given in Figure 15. Duchon's method with 


basis functions a? and a“ 


logd gave only a fair representa- 
tion of the bottom with regularly spaced model points. 
Additional model points did not improve the results so 
this technique was rejected. The basis function d logd, 
was introduced which gave good results (comparable to the 
multiquadric method) in one case and poor results in another. 
This was due to a dependence on the horizontal scale of the 
data. Independence of scaling for hydrographic survey 
modeling is very important since surveys are plotted at 
various scales. The method with basis function d logd is 
unacceptable for this reason. 

Shepard's formula gave best results in modified form 
with about six model points in each radius of influence. 
The inverse distance weighting function was better than 


the square of the inverse distance. The results were 
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considerably worse than those of Duchon's or Hardy's methods. 
Improvement could not be gained by increasing the number of 
model points. For these reasons, Shepard's method was 
unacceptable. 

Hardy's multiquadric analysis with inverse hyperboloids 
gave very poor results. For small Ô holes were produced 
at the model points and the depths between model points 
were not accurate. 

Of the methods tested, only multiquadric analysis with 
conic or sharply pointed hyperboloid kernels gave results 
which Indicated that further tests were warranted. Depic- 
tion of detail is improved by adding more model points with- 
out adversely affecting the model in other areas and the 


method is independent of linear scaling. 
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VI. FURTHER TESTING OF MULTIQUADRIC ANALYSIS 
WITH CONIC AND HYPERBOLOTD KERNELS 


The results from the previous section showed that 
multiquadric analysis with conic or sharply pointed hyper- 
boloid kernels was the only method tested that could meet 
the requirements of this application. Additional experi- 
mentation was done to determine the best procedures for 
selecting the model points and for joining models together 
at the boundaries. Tests were also run to determine how 
accurately the data sets could be represented with addi- 
tional model points while still saving significant storage 


space. 


A. SELECTION OF MODEL POINTS 


The selection of the data points to be used for the 
modeling is a critical process in the development of the 
multiquadric model. Three models for selection of the 
points were tested using the Auke Bay, Alaska data set. 
All point selection was done manually but consideration 


was given to the difficulty in automating the process, 


va Regular Spacing Selection 


In this method, data points from the survey were 


chosen at nearly even spacing without regard to depth, 
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bottom features, contour separation or any other factor. 

To avoid biasing, they were selected from a plot of record 
numbers rather than a plot of depths or depth contours. 
Additional points for more detail and accuracy were chosen 
for subsequent runs maintaining even spacing as much as 
possible without considering any factor except the hori- 
zontal distribution. The results of this procedure are 
presented in Table IX, The RMS differences were improved 
Significantly when the number of model points was increased 
from 53 to 110 but the maximum positive differences were 
not improved. Additional densification of the model points 
produced little improvement in either the RMS differences 


or the maximum differences. 


2. Iterative Selection 

In the iterative selection process, the results 
of one model were used to eliminate some model points and 
select additional ones to produce a better model. After 
developing the model with the first set of points, the 
comparison of survey data points with the model was analyzed. 
Additional model points were selected wherever single 
point comparisons showed the largest differences or in areas 
where several points showed relatively large differences 
of the same sign, Model points which had very small asso- 
ciated coefficients were eliminated. A small coefficient 
indicates that the associated basis function has little 
effect on the model since it remains near zero within the 
modeled area. The model was then computed with the new set 


BÉ points. 
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TABLE IX - Selection of Model Points with Even Spacing 


Number Maximum Maximum Number 
of model RMS positive negative of data 
points NPPR difference difference difference points 
55 O 2.02 Fel -5,64 1407 
Í I 200 9.45 -5.78 1407 
>> 2 1 91 a -6,39 1407 
52 10 190 SO 7,03 1407 
52: 15 aoe 8.68 -7.66 1407 
ÐI 20 195 3.68 -8,10 1407 
110 O 1537 9154 -4,25 1407 
110 Il 155 9,99 -4.71 1407 
110 2 150 10,094 -5.02 1407 
LO 5 126 10215 -5.46 1407 
TIO í 126 LORI -5.00 1407 
TLO 10 1226 10,28 -5.76 1407 
DO 15 120 10346 -5.97 1407 
144 li 1925 32368 -4.26 1407 
144 3 121 995 -4,66 1407 
144 2 1515 10307 -5.41 1407 
144 7 ryeal TORON -5.58 1407 
144 10 1.10 0200 -5.76 1407 
144 15 Tell TOO -6.01 1407 
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This procedure could be repeated until the desired 
accuracy was attained, the maximum number of model points 
to be used was reached, or the model accuracy no longer 
improved with further iterations, For this comparison of 
- selection methods, the process was repeated until the number 
of model points was approximately the same as the maximum 
number used in the test of the regular spacing selection 
method, 

Table X shows the results of these tests. In all 
tests, the best results were obtained when NPPR=0 (conic 
kernels). Both RMS and maximum differences improved sig- 
nificantly as the selection process was repeated. Two 
iterations yielded approximately the seme number of model 
points as the maximum used in the regular spacing selection 
method. 

Points related to features such as peaks, deeps 
or sharp changes in slope were chosen for the initial set 
of model points in these comparisons. Regular spaced points 
for the initial set were used in other tests with the 
iterative method. The results were good for both methods 
of initial selection. After a few iterations relatively 
few points from the initial set remained so the initial 
point selection method made little difference. 

A comparison of model point selection by regular 


Spacing and by iteration is shown in Figure 16. 
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TABLE X - Selection of Model Points by Iteration 


Number Maximum Maximun Number 

of model RMS positive negative of data 

points NPPR difference difference difference points 
82 O 1.48 4. 2 -5,52 1407 
82 1 1.56 2309 -5,50 1407 
32 5 2556 3.09 -9,43 1407 
107 O 1606 5.04 -3,29 1407 
107 Ji 1315 2701 -3, 32 1407 
LO] > 1,50 4,89 -4,52 1407 
152 O 0,72 1.91 -2.53 1407 
152 T O. 2.03 =2. 12 1407 
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Figure 16 - Comparison of Model Point Selection Methods 
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3. Complete Selection by Topographic Feature 


While using the iterative selection process, it 
was found that the additional points were selected where 
there were Significant changes of slope or where there 
were large areas without any model points. This led to 
an attempt to select all the model points in one step 
based on the following criteria. | 

e Select points at peaks, deeps, ridges and where 
slopes change significantly. 

e Select points to avoid leaving any large areas 
without model points as a result of the first 
criterion. 

A test was done by choosing 145 points to model 
only half of the Auke Bay data set. The RMS difference 
was 0.88 and the maximum difference was -3,22, These 
results show that one-shot selection is not nearly as 
good as the iterative method and probably not much better 


than the regular spacing method. 


4. Summary 


The iterative selection method gave by far the 
best statistical results. It also required the most com- 
puter time. It would be adaptable to complete automation 
Since the method for point selection has little subjectivity 


involved, 
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The regular spacing method would be easier to auto- 
mate but it doesn't give good results when detail is 
required. The method of complete selection by feature 
doesn't give significantly better results than does regu- 


lar spacing and the method would be difficult to automate. 


B. MODEL JUNCTIONS 


Agreement between two or more surveys which have a 
common boundary or cover a common area is a very important 
check on the quality of the surveys. Similarly, agreement 
between the models representing the surveys must be main- 
tained to avoid ambiguity. The problem is even more acute 
when several models are joined together to represent a 
Single survey. It would be desirable to represent an entire 
survey with a single model but in many cases this would be 
difficult due to the large systems of equations which would 
have to be solved to generate the model coefficients. 

Hardy (1971) suggested a simple method of junctioning 
where common points on the boundary were used in adjacent 
models. This would assure that the models were in agree- 
ment at these points and if ‘chosen at close intervals, the 
differences at intermediate points on the boundary would 
be relatively small. That method is not appropriate for 
this application since the data points are not along 
Straight lines which could be used as boundaries. Common 
points on irregular boundaries could be used but this would 


complicate model boundary definition and storage. 
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Two other more appropriate methods were r A 
for this application. Both use overlapping areas in the 
adjacent models rather than a common boundary. 

The first method is analogous to the method of using 
common points on a boundary because the model points within 
the overlapping area are required to be the same for both 

models. A line in the middle of the overlapping area is 


chosen to delineate the areas of model usage. See figure 17. 





- —— model A >+ 
hÜ —— model ES | 


k—use model 1—k— use model BE —~ 


x - model points for model A 
o =- model points for model B 
a ~ model points common to models A and 3. 


Figure 17 - Model Junctions by Overlap 


This method could produce small discontinuities at the 
boundary. 

The second method eliminates the discontinuity completely 
but requires more computation. The model points in the over- 


lapping area are not required to be common in both models. 
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In this area, a weighted sum of the values obtained from 
each model is used. The weight, w, is determined by the 
Hermite Polynomial 

w = 1 - 35° + 25” 
where s is the relative distance from the point of compu- 
tation to the overlap area boundary. The value of s varies 
from one at the outer model boundary to zero at the inner 


boundary of the overlap area (see Figure 18). 
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Figure 18 = Model Junction by Hermite Polynomial 
| 


The weight assigned to the model A value at the point of 
computation in Figure 18 is determined by using 


nad 5 
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The weight assigned to the model B value at the same point 


is determined by using 


a 
s = /D : (19) 
The sum of the resultant weights is always one. A plot of 


the Hermite polynomial is shown in Figure 19. 
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0 1 


Figure 19 - Hermite Polynomial 


Because the first derivative of this function is zero at 
s=0 and s=l, the transition from one model to an adjacent 
one will be smooth and continuous. 

The Auke Bay survey was divided into two overlapping 
models as pictured in Figures 17 and 18. Each was modeled 
Separately using the iterative method of model point selection, 
The two models were then joined by the two methods dis- 
cussed above. The results are presented in Table XI. 

Both junction methods showed improved results over the 
individual models since some of the largest errors were 
located near the outer boundaries of the models. The results 
with the polynomial method were only slightly better. Con- 
tour comparisons between the jumction methods showed little 
difference. The possible discontinuity at the boundary when 
not using the Hermite polynomial was not apparent in the 


contours, 
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These results have shown that transition from one model to 
another can be done smoothly. The method using the Hermite 
polynomial is only slightly better than the method using 
common points in an overlapping area. When joining models 
without common points, e.g. two different Surveys, the 
Hermite polynomial method will give a smooth non-ambiguous 


transition from one model to another, 


C. MODEL REFINEMENT 


Multiquadric analysis and the iterative method of model 
point selection were used to refine the models of all four 
data sets. Tables XII, XIII, XIV and XV give the statisti- 
cal results of these tests. Figures 20, 21, 22 and 23 show 
the effect of each iteration on the RMS differences, The 
number of model points is expressed as a percentage of the 
number of data points represented. This is a direct indi- 
cation of the storage savings attained by each model, All 
four figures show a similar trend. Initially, the results 
improve rapidly as the percentage of data points used in 
the model increases. The improvement then tends to level 
off and repeated iterations generate less improvement. 

Contour plots of the final models for each data set 
are shown in Figures 24, 25, 26 and 27. Comparison of 
these with Figures 7, 8, 9 and 10 show the agreement with 


contours of the original data. 


86 





For the Monterey Bay model, an RMS difference of 0,2 
fathoms and a maximum difference of -0.74 fathoms are good 
considering the irregularity of the bottom. These results 
were obtained using 17.2% of the data points as model 
points. 

For the Morro Bay model, an RMS difference of 0,55 
feet and a maximum difference of -1.58 feet are good consikr- 
ing that the Neots range from 16 to 82 feet. This repre- 
sentation was made using only 15.7% of the data points. 

The allowable error specification (Section I.B) is one foot 
in the shallow end of this range and three feet in the deep 
end. 

For the Auke Bay data set, an RMS difference of 0,530 
fathoms and a maximum difference of -0.95 fathoms using 
20.6% of the data points are good considering the steep 
slopes in the area. A horizontal positioning error of a 
few meters (within tolerance for the survey scale) could 
create several fathoms difference in many of the recorded 
depths. 

Even though the range of depths in the Gulf Coast data 
set is small, and RMS difference of 0.30 feet and a maximum 
difference of -1.16 feet using 9.9% of the data points 
could be considered good. The allowable error (Section I.B) 
for these depths is one foot. There are places in the data 
set where crossline and adjacent soundings from multiple 


vessels disagree by as much as two feet. The model 
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representation tends to smooth out such discrepancies and 
this smoothing appears as relatively large differences in 
the statistics. 

The Gulf Coast model gave the only case where a large 
Ò (WPPR=50) gave much better results than smaller values. 
Only nine model points were used in that case. When more 
model points were added a small NPPR was required for good 


results. 
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Figure 20 - Monterey Bay Model Results 
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Figure 21 - Morro Bay Model Results 
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Figure 22 - Auke Bay Model Results 
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Figure 23 - Gulf Coast Model Results 
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Figure 24 = 226 Point Monterey Bay Model Contours 
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Pigure 27 - 165 Point Gulf Coast Model Contours 
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VII. CONCLUSIONS 


Of the methods examined, only Hardy's multiquadric 
modeling technique with conic kernels was found to be suit- 
able for modeling hydrographic survey data. Polynomial, 
double Fourier series and finite element methods were 
rejected for reasons found in the literature. Duchon's 
method, Shepard's formula, and multiquadric analysis with 
inverse hyperboloid kernels were rejected as a result of 
tests presented herein, 

Selection of model points for the multiquadric method 
is best performed by iteration. An initial set of model 
points may be chosen either by regular spacing or by bottom 
feature. Comparisons of the initial model with the survey 
data are used to select additional model points where there 
are large errors. Points with the smallest coefficients 
are eliminated from the set. A new model is computed and 
the process is repeated. This procedure can be repeated 
until the desired accuracy is reached or until additional 
iterations fail to produce increased accuracy. This metnod 
is considerably better than selecting regular spaced points 
and densifying them for more accuracy or selecting the maxi- 
mum number of points at one time by examining the bottom 
features. The iterative method of model point selection is 
adaptable to automation since there is relatively little 


subjectivity involved. 
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Adjacent survey models or partial survey models can 
be joined to produce a continuous and unambiguous repre- 
sentation of the bottom in the junction area. The best 
method is to use a weighted sum of the model values in an 
overlapping area. The Hermite polynomial function should 
be used to produce the weights. If the two models contain 
common points in the overlapping area, a good junction can 
be made without computing a weighted sum. This method 
could produce a slight discontinuity along the boundary 
but it was not apparent in the test runs for this project. 

The multiquadric technique with iterative selection 
of model points and Hermite polynomial junctions produced 
good models of four data sets. Approximately 20% of the 
survey data points were required for a model of Auke Bay, 
Alaska, where there is tremendous bottom irregularity. 
Only 10% of the data points were required for a model of 
the Gulf Coast where the bottom shows little variation. 
Approximately 15% and 17% were required for the Morro Bay 
and Monterey Bay models where there is more irregularity 


and moderately sloping bottoms. 
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